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Abstract Both community ecology and conservation biology seek further understand- 
ing of factors governing the advance of an invasive species. We model biological invasion 
as an individual-based, stochastic process on a two-dimensional landscape. An ecolog- 
I— { , ically superior invader and a resident species compete for space preemptively. Our 

general model includes the basic contact process and a variant of the Eden model as 
Q\ \ special cases. We employ the concept of a "roughened" front to quantify effects of dis- 

creteness and stochasticity on invasion; we emphasize the probability distribution of 
the front-runner's relative position. That is, we analyze the location of the most ad- 
vanced invader as the extreme deviation about the front's mean position. We find that 
a class of models with different assumptions about neighborhood interactions exhibit 



O . universal characteristics. That is, key features of the invasion dynamics span a class 



of models, independently of locally detailed demographic rules. Our results integrate 
theories of invasive spatial growth and generate novel hypotheses linking habitat or 
O^i landscape size (length of the invading front) to invasion velocity, and to the relative 

position of the most advanced invader. 
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1 Introduction 



C^ ' Invasive species threaten biodiversity in many ecological communities (|Ruesink et al. 19951 

C^ , |Rosenzweig 200l| , and impose significant costs on agriculture and public health ifPimentel et al. 2000)l . 

Spatial advance is the most obvious phase of invasion ([Shigesada and Kawasaski 1997| 
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[Hastings et al. 2005[ ), and perhaps the most important (jLockwood et al. 2007|l . In- 
deed, the dynamics of advancing fronts challenges our capacity to predict how in- 
vasive species, emerging infections, and evolutionary adaptations spread across land- 
scapes (|Andow et al. 19901 |0'Malley et al. 2006a[ ). Our study emphasizes certain es- 
sential features, signature characteristics, of advancing fronts that let us integrate some 
common models of stochastic spatial growth. 

We model preemptive competition between an invader and a resident species as 
a stochastic process in a two-dimensional environment. As the invading competitor 
advances, the model exhibits variability of the invader's location about the front's 
mean position. Generation of variability along the front separating invader and resident 
species is called interface roughening (Krug and Meakin 1990 ; Barabasi and Stanley 1995| 
|Halpin-Healy and Zhang 1995[). In addition to invasion processes in population biol- 



ogy (O'Malley et al. 2006b I, roughened fronts emerge commonly in other systems, in- 
cluding the growth of cancerous tumors (|Bru et al. 2003p . task-completion landscapes 
in in massively parallel computer networks (|Korniss ct al. 2000; K orniss et al. 2003|) . 
and surface growth in materials science ( [Barabasi and Stanley 1995[ ). Models built on 
quite different assumptions can roughen in the same way; features dependent on rough- 
ening will then behave similarly across models. 

We recently analyzed the time course of roughening in our two-species model for 
competitive invasion ( [O'Malley et al. 2006b[ |. In this paper we focus on properties of 
the invader's extreme advance, the front-runner's location (lEllner et al. 1998)) . We in- 
vestigate effects of stochasticity, discreteness and habitat size on the invasion process 
(Habitat, or landscape, size refers to the linear extent of the system transverse to 
the direction of invasive advance). For particular parameter values and initial condi- 
tions, our model for spatial invasion becomes equivalent to the contact process (CP) 
([Harris 19741 IDurrett and Levin 1994allHinrichse n 2000; Obor ny et al. 2005[ ). Another 
parameter set produces a variant of the Eden model (,Eden 196LIJullien and Botet 1985al 
IJullien and Botet 1985b] IKawasaki et al. 2006[l : both CP and the Eden model have 
served as archetypes for modelling invasive spread. We aim to extract features com- 
mon to a family of models (hence, characteristics of a universality class). To do so, we 
first investigate the scaling behavior of our model's interface width (variability of the in- 
vader's advance about its mean position), both as a function of time and of habitat size. 
We then address the shape of the probability density of the front-runner's position. Our 
analyses reveal consistent universalities. Model details, of course, affect non-universal 
coefficients (or prefactors), but universal features lead to predictions spanning an im- 
portant class of stochastic growth models. Furthermore, although the front's velocity 
depends on details of a model's local birth and death rates, the approach to the front's 
asymptotic velocity - again, as a function of time and of habitat size - is universal, and 
these dependencies offer a useful tool for identifying ecologically salient properties of a 
given invasive front. 



1.1 Advancing fronts: background comment 

Spatial expansion characterizes biological invasion from local to geographic scales 
( [Minogue and Fry 1983[lvan den Bosch et al. 1992l[Sh igesada and Kawasaski 1997p . Se- 
lecting a framework for modelling spatially detailed, invasive growth has been posed 
as a choice between reaction-diffusion equations for continuous densities or stochas- 
tic discrete (individual-based) models ([Durrett and Levin 1994bl) . The former include 



deterministic reaction-diffusion systems, which often permit approximation of the in- 
vader's speed as the asymptotic velocity of a travelling wave ( |Dwyer and Elkinton 1995] 
ICaraco et al. 20021 [Murray 2003[ ). Significant generalizations of the basic determin- 
istic theory vary schedules of reproduction or the distribution of dispersal distance 
(jCantrell and Cosner T99T]|Shigesada et al. i995|IKot et al. 1996IINeubert and Caswell 2000] 
|Dwyer and Morris 2006[ ) and treat advancing fronts as a travelling waves with con- 
stant or increasing velocity. However, travelling waves invoke infinitesimal popula- 
tion densities (van Baalen and Rand 1998[l . The front can be "pulled" by growth and 
spread of the invader at locations where its density is near ( [Snyder 2003[ ); the model 
may consequently overlook realistic features of the dynamics of rarity ([Lewis 20001 
[Clark ct al. 2003 , Antonovics et al. 2006 ) . That is, deterministic reaction-diffusion equa- 
tions are not capable of capturing ecologically relevant effects of spatially correlated 
variability along the types of fronts we study. 

Nevertheless, stochastic partial differential equations (or Langevin equations) and 
equivalent field-theoretical descriptions have, for decades, been employed successfully as 
coarse-grained (i.e., "mesoscopic," continuous-density) representations of both interact- 
ing particle systems in physics and chemistry (IKorniss 19971 iPechenik and Levine 1999[l . 
and individual-based models in biology ([Moro 20011 lEscudero et al. 2004]| . In princi- 
ple, for every individual-based model, one could, in the continuous-time limit, con- 
struct a dynamically faithful coarse-grained representation using a stochastic partial 
differential equation ( [van Kampen 1981[ IGardiner 19851 ISchmittmann and Zia 19951 
iHinrichsen 200Cl|. Deriving the appropriate noise terms presents a challenge, and the 
resulting model, with some exceptions, will not likely yield an analytically tractable 
solution. One then often resorts, as in the case of individual-based models, to numerical 
investigation. The important point, for questions we address, is that properly developed 
Langevin equations ( [van Kampen 1974[ IDoi 1976] IPeliti 1985[) and spatially explicit, 
individual-based models should exhibit identical scaling behavior when critical, large- 
scale levels of variability govern properties of the invasion dynamics. We therefore 
choose the more familiar individual- based approach ( [DeAngelis and Gross 1992[ ). 

Individual-based constructions capture effects of nonlinearity and stochasticity gov- 
erning invader dynamics at introduction and at an advancing front ([Wilson et al. 19931 
I Wilson 19981 Irhomson and Eli nor 2003). For both individual- based and stochastic par- 
tial differential equation models, if diffusion is limited, the front is usually "pushed" at a 
velocity less than that of the corresponding deterministic diffusion model ([Moro 2003p . 
Furthermore, fronts in these stochastic models roughen; that is, they fiuctuate strongly 
([Kardar et al. 19 86 ) . Given an initially fiat interface between species, the invader quickly 
advances different distances along the length of the front. As the front roughens, the size 
of random fiuctuations along the front become spatially correlated ([Racz and Galfi 19881 
Iben Avraham 19981 [Majumdar and Comtet 2004[ ), and this correlated variation ex- 
erts significant, time-dependent effects on the front's velocity ( [Krug and Meakin 1990[ 
Ivan Saarloos 2003)) . Roughening also implies that the time-dependent position of the 
invader's furthest advance, the "front-runner," will be influenced by the same spatial 
correlations. One-dimensional models cannot, of course, generate a roughened front. 

Our analysis of the front runner's position involves the extreme value among 
strongly correlated random variables; Fisher and Tippett (1928) and Gumbel (1958) 
first derived the distribution of the extreme value among independent random vari- 
ables. For roughened interfaces, results for independent (or weakly-correlated) random 
variables break down because of the strongly correlated levels of variability along the 
front. Recently Majumdar and Comtet (2004, 2005) analytically obtained the scaling 



behavior of the distribution of the extreme advance of the Kardar-Parisi-Zhang inter- 
face (|Kardar et al. 1986|) . The interface-set that emerges in invasion models studied 
here belongs precisely to this universality class; hence one immediately obtains the 
shape of the probability distribution of the extreme advance (the front-runner) in the 
steady state. 



2 Model 

Invaders often must overcome competition (|Elton 19581 IParker and Reichard 19981 
|Hoopes and Hall 2002 1 ISimberloff et al. 2002p to advance spatially. We assume pre- 
emptive competition; neither species can colonize an occupied site until the occu- 
pant's mortality opens that site (IComins and Noble 19851 [Connolly and Muko 2003] 
ITainaka et al. 20041 lYurkonis and Meiners 2004[l . Since we address frontal advance, we 
ignore introduction from beyond the habitat; a species may occupy new sites only 
through local propagation when one or more of its nearest-neighboring sites is empty. 
Although some invasions may be driven by long-distance dispersal (|Ferrandino 19961 
ILewis 19971 IClark et al. 200l|) . our assumption of limited dispersal has clear empiri- 
cal justification. D'Antonio (1993) reports that 98 per cent of an invasive succulent's 
{Carpobrotus edulis) population growth over four years resulted from spatially clustered 
clonal propagation. Invasive Argentine ants {Linepithema humile) advance spatially by 
locally budding new from existing colonies ( [Holway 1998[ ). For many invasive plants, 
and some invading animals, phalanx-like advance should prove realistic within any 
single habitat. Importantly, we note that adding explicit finite-range diffusion, with a 
small rate, to local colonization does not change the universal aspects of the quantities 
studied here (|Moro 20011) . 



2.1 Model dynamics and simulation 

On an Lx x Ly lattice, a site represents the minimal resources necessary to maintain a 
single individual. The local occupation number at site x is nj(x) = 0, 1 with z = 1, 2, 
referring to the resident and invader species, respectively. During a single time unit, 
one Monte Carlo step per site [MCSS], LxLy sites are chosen randomly for updating. 
An empty site may be occupied by species i through propagation from a neighboring 
site at rate airii{x.), where a^ is the individual-level propagation rate for species i, and 
rji{:>i) — (l/5)^x'enn(x)^j(^') i^ the density of species i in the neighborhood around 
site x; nn(x) is the set of nearest neighbors of site x, and 5 is the number of sites in 
that neighborhood {S = |nn(x)|). Unless noted otherwise, the results presented here 
are for 5 = 4, but in a few cases we set (5 = 8 (Moore neighborhood) and 5 = 12 (Moore 
neighborhood plus 4 sites) for comparison. An occupied site opens through mortality of 
the individual. We assume density-independent mortality (|Cain et al. 1995")) and assign 
each species the same mortality rate ^. Summarizing transition rules for an arbitrary 
site X, we have 

-^ ' 1, -A ^ 2, 1-^0, 2-^0, (1) 

where 0, 1, 2 indicates whether a site is open, resident-occupied, or invader-occupied, 
respectively. Table [l] defines the symbols and notation we use. 



Table 1 Definitions of model variables and parameters 





Symbols 


Definitions 




Lx, Ly{=L) 


Lattice size 




X 


Location of lattice site 




ni(x) 


Occupation number for residents at site x; ni (x) = 0,1 




n2(x) 


Occupation number for invaders at site x; ?12 (x) = 0,1 




nn{x) 


Set of nearest neighbors around site x 




<5 


size of neighborhood around site x (5 = |nn(x)|) 




»?i(x) 


Density of species i on nn(x) 




ai 


Individual rate of propagule production, species i 




H 


Mortality rate 




ac{fJ,) 


Minimal propagation rate for persistence 




hy{t) 


Rightmost invader in row y at time t 




h(t) 


Mean of hy{t) 




fT'maxVt) 


Rightmost invader at time t 


Amaxit) = hmax{t) - h{t) 


Distance from front-runner to mean of front 




V* 


Asymptotic velocity of invasive front 




{w^) 


Mean squared interface width 




m 


Correlation length along front 




a 


Roughness exponent 




P 


Growth exponent 




z 


Dynamic exponent 




P^* 


Equilibrium single-species density 




N" 


Number of invaders in a row through the width w 



We assume that interspecific competition drives the dynamics (i.e., each species 
persists absent competition), and that the invader has a reproductive advantage. There- 
fore, we restrict attention to the ac{iJ^) < ai < a2 regime, where ac(M) is the critical 
propagation rate below which either species, in the other's absence, grows too slowly 
to avoid extinction ( [Oborny et al. 2005[|Q'Malley et al. 2006al ). 

We impose periodic boundary conditions along the y-direction of the Lx x Ly lat- 
tice. The initial condition is a flat linear front (straight vertical line), i.e., the invader 
completely occupies a few vertical columns at the left edge of the lattice, and all re- 
maining sites are occupied by the resident. Invasive advance, therefore, proceeds in 
the a;-direction. As a simulation begins, mortality quickly reduces the density on each 
side of the front to the respective "quasi-equilibrium" value, where a species' propa- 
gation balances its mortality. As a simulation continues, we track the location of the 
invading front by defining the edge as the location of the right-most individual of 
the invading species, hy{t), for each row y (Fig. [!}. We shall also refer to this quan- 
tity as the local "height" [borrowing terminology from the non-equilibrium surface 
and interface-growth literature ( [Barabasi and Stanley 1995[ )]. We record the average 
position h{t) — (i/Ly)'^ hy{t) for each time step. We estimate velocity once h{t) 
approaches linear increase with time. We ran each simulation until the front reached 
the end of the system. Longitudinal system size Lx has no particular impact on system 
behavior. But the transverse system size Ly plays a fundamental role in controlling the 
large-scale properties of the emerging fronts and the interface region. Hence the size of 
the habitat (specifically, the length of the front) will exert an important influence on 
the dynamics of ecological invasion. 

If the resident is absent initially, model assumptions preclude its introduction. Elim- 
inating the resident species (i.e., ni(x) = for all x) reduces the model to the basic 
contact process for invaders. Furthermore, for the specific choice of /x = (i.e., no mor- 
tality), the CP itself reduces to a variant of the Eden model (|Jullien and Botet 1985a]) . 



t=]00 



t=1000 



t=4000 
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02=0.80, fi=0.10. Time has units of Monte Carlo steps per site (MCSS). White cells represent 
empty sites, while gray and black correspond to sites occupied by residents and invaders, 
respectively. 



To emphasize the strength of the universal aspects of these stochastic growth processes, 
we also present some results on the front's extreme advance for these special parameter 
choices. 

An advancing front "roughens" (Fig.[l]). That is, the typical size of deviations about 
the mean front position increases until the distribution of the mean squared fluctua- 
tion equilibrates statistically. The resulting steady-state depends on habitat size Ly 
(|Barabasi and Stanley 1995| |Halpin-Healy and Zhang 1995[ ) . Figure [2] shows a typical 
configuration with "width" w about the front's average position h{t). Similar "wavy" 
front configurations can be observed in results plotted in recent front-velocity stud- 
ies of the basic CP (|Ellner et al. 1998)) and of the Eden model (|Kawasaki et al. 2006^ . 
This "waviness" is at the heart of kinetic roughening, the focus of our investigation. 
We emphasize that as habitat size Ly increases, the amplitude of the interface fiuc- 
tuations (deviations about the mean) in the steady state increases monotonically. In 
turn, one can employ universality arguments and recent results to find the shape of 
the distribution of the front-runner's relative position ( |Majumdar and Comtet 2005[ ) 
and important habitat-size dependent corrections to the velocity of invasive advance 
dKrug and Meakin 19"90l ). 



3 Roughening of the front 

Next we apply a framework for dynamic scaling of roughened interfaces to advance a 
new perspective on invasion ecology. To analyze roughening we define the "width" of 
the invading front via: 



Tuu2 



w'{Ly,t) = — \[hy{t)-h(t)] 



(2) 




Fig. 2 Width (ui) and the extreme advance {Amax) relative to the mean front position (h) 
in a rough front. For illustration, correlation length ^ is also indicated. 



Of course, w {Ly,t) is itself a fluctuating quantity, so in simulations we estimate an 
ensemble average over different runs, or take an average over a time series in steady 
state, to obtain {w {Ly,t)). Hereafter, we let L = Ly, for convenience, and write 
{w {L,t)). Without creating confusion, we shall use the notation w — y (w^(L,t)), to 
refer to the typical size of the interface region. 

w {L,t) in Eq. ([2]) is a function of the local stochastic invasion heights hy{t), and 
the hy (t) axe not independent random variables. Starting with a flat interface, a single 
correlation length ^(t) develops along the front (Fig. [5|. This correlation length grows 
as a power-law function of time ^(f) ~ f ' ^ ( [Barabasi and Stanley 1995[ ), but once ^(i) 
spans the length L of the system, at a crossover time on the order of tx ^ L^ , rough- 
ening has reached its stationary state. The heights of invasive advance along the front 
have become and remain dependent random variables, complicating analysis of the 
invader's maximal advance ( |Raychaudhuri et al. 200"T] |Majumdar and Comtet 2004[ ). 
We let (w [L, oo)) represent the mean squared deviation when its distribution has 
equilibrated; below we focus our front-runner analysis on these steady-state fluctua- 
tions. 

Advancing fronts with different locally structured propagation and mortality pro- 
cesses may still exhibit the same dependence of roughening on time, and the equi- 
librium variability along the front may exhibit the same dependence on habitat size. 
Such fronts belong to the same "universality class"; universality offers powerful sim- 
plification and generalization ( |Cardy f996j Ferreira and Alves 2006 1 ). O'Malley et al. 
(2006b) found that our model's roughening behavior belongs to the KPZ universality 
class (for Kardar-Parisi-Zhang; see Kardar et al. 1986). We summarize those results, 
since quantifying roughening allows us to correct the asymptotic velocity of invasion, 
and to write the probability density of the front runner's relative position. Appendix B 
provides some further details about the KPZ equation, a continuous stochastic differ- 



ential equation, that captures interface dynamics for a large class of discrete, stochastic 
individual-based systems. 

The profound nature of roughened fronts (or self-afEne interfaces) lies in the inti- 
mate connection between their temporal and spatial scaling behavior (S arabasi and Stanley 1995| 
|Halpin-Healy and Zhang 1995[ ). For a system with transverse system size L, for early 
times, the width (or spread) of the interface [Fias. [Tl2] increases in a power-law fashion 
{w {L,t))'^t P . Then, after a system-size dependent cross-over time tx^L^, the front 
reaches its steady-state: interface fluctuations about the steadily advancing mean front 
become stationary, but their magnitude will depend on the transverse system size. In 
particular, the interface width "saturates" in time, and its stationary value will scale 
with the system size as {w (Z/,oo))~L ". Thus, both the time to reach steady state 
and the size of the steady-state width depends on the transverse system size L in a 
power-law fashion. Put simply, if the transverse system size is increased, both the time 
to reach steady state and the width of the interfacial region in the steady state in- 
crease. The exponents q>0, /3>0, and z>0 are referred to as the roughness, growth, 
and dynamic exponent, respectively. (The roughness exponent a should not to be con- 
fused with the local propagation rates of the two species.) These fundamental scaling 
properties of roughening can be summarized as 

/ 2/r ,N\ fi^*^ f0rt<fx , rZ /n^ 

(»(i,i)>-|^2.fo^,>>,^ ' *x~^ (3) 

and are illustrated schematically in Fig. [31 In general, these exponents obey a scaling 
relationship a — j3z ( [Barabasi and Stanley 1995[ ). 

Models whose kinetic roughening can be described by the same set of exponents arc 
said to belong to the same universality class ( Barabasi and Stanley 1995|[Halpin-Hcaly and Zhang 1995[ ) 
Fronts generated by models with different local rules often belong to the same universal- 
ity class (for example, fronts in the Eden model, the contact process, and the two-species 
model of this paper). Importantly, models formulated to address completely different 
questions, ranging from surface physics and computer science (|Korniss et al. 2003|) to 
ecology, can belong to the same universality class. 

The front in our two-species individual-based model exhibits scaling consistent with 
Expression [31 the exponents are /3=l/3, a=l/2, and z=3/2 ( |0'Malley et al. 2006b"| ). 
Consequently, the front belongs to the so called KPZ-interface universality class (in 
one transverse dimension) (Fig. [J). Slight deviations from the theoretical exponents 
result, in part, from corrections to scaling and "intrinsic width" effects for limited 
system sizes (jMoro 2001l|Blythe and Evans 2001[ ). Furthermore, the strong correlation 
in the steady-state time series of w (t) (the auto-correlation time also diverges with 
the system size as L^) increases sampling error for larger system sizes. For the KPZ 
universality class, the analytical result by Foltin et al. (1994) specifies the universal 
shape of the distribution of the steady-state width (see Appendix B.2). The results 
can be used to identify the universality class of a model with a rough front; doing so 
provides further support that our model belongs to the KPZ class (Appendix B.2). 



3.1 Roughening and the asymptotic invasion velocity 

Estimating the speed at which an invasive species advances remains a focus of ecological 
theory l|Mollison and Levin 1995llEllner et al. 1998llLewis 2000IIKawasaki et al. 2006p . 
Finding an analytical expression for a stochastic front's asymptotic velocity v* presents 
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Fig. 3 Schematic plot illustrating tiie scaling features of the width of rough fronts, Eq. ((Sjl- 
(a) Interface width vs time on log-log scales, for different system sizes Li < L2 < L3. For early 
times, the width of the interface increases as {w'^{L,t))^t'^P . After a system-size dependent 
cross-over time ix {L)^L^ , the interface fluctuations become stationary, and the width reaches 
its system-size dependent steady-state value {w^ {L, oo))'^L^°' . The straight solid line corre- 
sponds to the power- law increase for early times before the crossover (growth phase). For larger 
transverse system size L, it takes longer to reach steady state, and the width in the steady 
state is larger, (b) Scaled plot {w^{L,t))/L'^°' vs t/L^, yielding full data collapse for different 
system sizes, the signature of dynamic scaling jBarabasi and Stanley 1995| l. The straight solid 
line indicates the power-law behavior x'^" of the scaled variable x=t/L^ for x <^ 1. 



a challenge (|Pechenik and Levine 1999)) . particularly for two-dimensional fronts. Ana- 
lytically tractable, pulled fronts emerge naturally in the deterministic reaction-diffusion 
limit of our model ([Fisher 1937IIHolmes et al. 1994l[Diryer and Elkinton 1995| [Murray 2003] 
Ivan Saarloos 2003| (see Appendix A) . But their velocity typically overestimates that of 
the actual individual-based model in the diffusion-limited regime l|Moro 2001IIMoro 20031 
[Doering et al. 2003[ ). Our model, with only local dispersal and no "explicit" diffusion, 
lies precisely in that regime: since interspecific mixing is insufficient about the in- 
terface, and internal fluctuations are important, deterministic equations break down 
([Antonovics et al. 200"6] |0'Malley et al. 2006"bt . 

Deterministic reaction-diffusion models neglect discreteness of individuals, the fun- 
damental source of endogenous fluctuations (|Escudero et al. 2004|) . Consequently, de- 
terministic reaction-diffusion theory neglects the substantial, correlated random vari- 
ation in the extent of the invader's advance along the front. To demonstrate the sig- 
niflcance of this difference, we compared results from Monte Carlo simulation of the 
individual-based model defined by the transition rates in Eq. ([T| to the velocity of 
travelling waves (pulled fronts) arising in our model's deterministic reaction-diffusion 
limit, Vp = {^J./Q.l)^/ 0-2(02 — cti) (see Appendix A). Asymptotic invasion velocity Vp 
increases as either the common mortality rate fi or competitive asymmetry, [02 — Qi), 
increases; a strong resident competitor can slow an invasive species' advance under a va- 
riety of assumed dynamics (|Hosono 19981 IGandhi et al. 19991 [Weinberger et al. 2002| 
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Fig. 4 (a) Average width as a function of time (log-log scales) for various system sizes, 
averaged over 20 independent realizations with oi=0.70, 0:2=0.80, and ^=0.10. Dashed line 
corresponds to the one-dimensional KPZ power law with the exponent 2/3 = 2/3. Inset shows 



the scaled plot, {w^{L,t))/L^ 



t/L'' 



using the KPZ exponents, (b) Steady-state width, 
w = a/(ui2(L)), and extreme fluctuations, {Amax(L)), (averaged over time series in the steady- 
state) as a function of system size, L{= Ly), for two sets of parameters: Set 1 (q:i=0.50, 
02=0.70, and /i=0.20); and Set 2 (oi=0.70, 02=0.80, and /i=0.10). Dashed lines indicate the 
best-fit power laws. The solid line corresponds to the the theoretical asymptotic scaling in the 
KPZ universality class with roughness exponent a=l/2. 



|0'Malley et al. 2009[ ). For further comparison with our stochastic model, note that the 
traveUing wave approaches its steady-state value as v{t) = Vp — 0{l/t); that is, the 
travelling wave's temporal correction decays as t~ HHolmes et al. 1994p . 

Figure [Sfa) plots asymptotic velocities from simulation of our individual-based 
model (for fixed 0:2) as a function of competitive asymmetry (02 — Ql)- For pa- 
rameter values chosen, the invader advances more slowly than our model's deter- 
ministic reaction-diffusion limit (i.e., slower than Vp) [Fig. [SJb)]. It is important to 
note that the pulled front approximation Vp can fully reproduce the numerical solu- 
tions of the deterministic nonlinear reaction-diffusion system [Appendix A, Eqs. (|A.1|I 
and (|X2)l ] ( |0'Malley et al. 2006b| |Q'Malley et al. 2009| ). What it fails to account for 
are the effects of stochasticity. In our discrete, stochastic model the front becomes 
pushed (Ivan Saarloos 2003ll : invasion velocity depends on the full frontal region, with 
its nonlinear interactions, rather than depending on the leading edge only (a pulled 
front). For comparison, we also plot an approximate expression for the velocity, sug- 
gested by Moro (2003), vjji^ = ^(02 /«l ~ 1), applicable to strongly diffusion-limited 
systems (see next section), where the effective diffusion coefficient is negligible com- 
pared to the effective reaction rates. As Fig.[5l[b) shows, except for large (0:2 — Ol), this 
approximation underestimates the front velocity of our individual-based model. This 
may be expected by noting that the effective reaction and diffusion coefficients are of 
the same order of magnitude (see Appendix A) . 

The functional expression for a model's asymptotic velocity cannot be universal, 
since it must depend on details of the local dynamics. However, Krug and Meakin (1990) 
showed that temporal and finite-size {Ly — L) corrections to the asymptotic velocity 
are universal within a given class (see Appendix B.3); a general ecological prediction 
consequently emerges. Frontal velocity initially increases until reaching steady state. 
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Fig. 5 (a) Asymptotic front velocities from Monte Carlo simulations of the individual-based 
stochastic model for fixed 02=0.70, as a function of the difference between propagation rates, 
0:2— ai, with Ly = 100, for three values of fi. (b) Simulation results from (a) for fi = 0.20 com- 
pared with the results from numerically iterating the deterministic nonlinear reaction-diffusion 
equations. Appendix A, Eq. i TOl (O'Mal ley et al. 20 06b O'M alley et al. 2009 ) and with two 
approximate analytic limits: pulled-front approximation [solid curve. Eg. | |A.4|| ] and diffusion- 
limited front [dashed curve, IIMoro 2003| |]. The solid straight line is the best-fit effective power 
law in the region where the difference between propagation rates is small, corresponding to 
I) ~ (02 - ai)* with = 0.66 ± 0.04. 
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Prior to the steady state, the difference between the current and asymptotic velocities 
(the temporal correction) scales similarly with time for every member of a particular 
universality class. More importantly, the second correction recognizes that steady-state 
velocity depends on habitat size; the habitat-size correction is proportional to L~ 
[see Appendix B.3, Eq. (|B.6p ]. That is, a given invasive species' front will advance 
faster as linear habitat size increases. More generally, the finite-size correction to the 
asymptotic velocity scales similarly with habitat size for each stochastic model where 
invasive advance exhibits a KPZ interface. 



3.2 Invader spatial profiles and diffusion limitation 



Next we focus on invader density within the width of the advancing front, to compare 
stochastic and deterministic invasion models in a general context. We constructed 
density profiles for the invader along the horizontal direction (averaged over rows y, 
for a given x): 



P2{x) = —^n2{x,y) 



(4) 



The density profile should spread horizontally in response to the front's roughening 
(as a function of time during growth phase, and as a function of habitat size in steady 
state). Figure [6l^ a) shows the invader's steady-state density profiles as a function of 
distance from the front's mean position. Profiles are averaged over 10 time steps, 
and plotted for a series of habitat sizes L. Behind the front, where the resident has 
been excluded, invader density p2{x) rests at its single-species equilibrium p2* ■ Invader 
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(x - h)/C 

Fig. 6 (a) Steady-state spatial density profiles of the invader species [averaged over all rows, 
Eq. (|4]|] for several habitat sizes, for parameter "Set 1" (oi=0.50, a2=0.70, /i=0.20), averaged 
over 10* time steps, (b) Same data as in (a), but scaled according to the system-size dependent 
width of the front {w ^ L°' with a = 1/2). 



density declines through the width, and only the resident occupies locations well right 
of the front. 

The steady-state density profiles show that the span of the frontal region, within 
which P2* > P2{^) > 0, increases progressively with habitat size L [Fig. [6l^a)]. Since 
steady-state roughening scales with habitat size according to {w {L, a)) ~ L ", rescal- 
ing the abscissa by L" {— L ' ) results in the different density profiles falling on the 
same plot [Fig.|6jb)], an articulate data collapse indicating the fundamental, underlying 
dependence (cf. Antonovics et al. 2006). 

We stress that expansion of density profiles with increasing habitat size is not a 
result of greater mixing of the two species within the span of the front. Rather, density 
profiles result from emerging spatial correlations between fluctuations in the invader's 
advance along the length of the front, a collective result of the neighborhood-level 
stochastic dynamics. Indeed, spatial configurations in Figs. \T\ and [2] reveal that few 
single rows, if any, suggest an invader-density gradient. In fact, the density profile in 
most rows resembles a "shock- wave" . 

Moro (2003) analyzes transitions from pushed-front to pulled-front velocities in 
discrete, stochastic models. Our simulations fall into Moro's diffusion-limited regime, 
where strong spatial correlations produced by local propagation dominate effects of 
mixing along the front. Our model's assumptions (in particular, no explicit diffusion) 
and parameter values imply that the mean-field assumption of strong mixing at the 
between-species interface cannot hold. Consequently, our competitive invasion pro- 
duces, in any given row, an invader-density gradient that spans only a narrow region; 
essentially, the invading species advances as a Shockwave front. Moro's (2003) analy- 
sis identifies the significance of A''* , the (approximate) number of invaders occupying a 
row of the interface between the two regions of single-species equilibrium densities. The 
front's velocity asymptotically approaches the reaction-diffusion (pulled-front) velocity 
Vp only as A^* — >■ oo. For a mortality rate oi fi = 0.1, our simulations have 1 < N* < 
2, so that the impact of discreteness and stochasticity remains strong. 
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4 The front-runner: roughening and scaling of extremes 

Finally, we turn to the problem of the invading species' "front-runner" (|Ellner et al. 19981 
IClark et al. 200D IThomson and Ellner 2003|) . and scaling of the distribution of the 
maximal invasive incursion. The front-runner's position presents an extreme value prob- 
lem, complicated by the probabilistic, strongly correlated dependence of the invader's 
advance at different locations along the front ( Raychaudhuri et al. 2001 Maju mdar and Comtet 2004[ ). 
Hence, traditional extreme- value statistics (Fisher 1928: Gumbel 1928 Berman 19641 
IGalambos 19871 IGalambos et al. 1994[l are not applicable. 

We restrict attention to height fluctuations occurring after roughening has equili- 
brated {t » tx^L^)- At each time step we identify the farthest-advanced location 
of the invading species, hmax(t) ~ inaxy{hy{t)}. Measuring the extreme advance rel- 
ative to the front's mean position, h{L,t), we obtain the maximal relative height at 
time t, Amax{t) — hmax{t) — h{L,t). We sample Amax{t) sufficiently to construct his- 
tograms to estimate the probability density of the extreme fluctuations P{Amax,L), 
for different propagation/mortality rates, and for different habitat sizes L. 

Results presented above led us to conclude that the front in our two-species invasion 
model belongs to the KPZ class. Hence, we can proceed to characterize properties of 
invasive extremes. Given steady-state roughening of the invading species' front, the 
probability density of the extremes can be written as 

P{Amax,L) ^ -— r<P{Amax/{Amax)) , (5) 

{•^max) 

where <f'(it) is the probability density of the scaled variable u = Amax/{Arnax) and 
is independent of habitat-size L. !Z'(w) is a scaled Airy density for all interfaces belonging 
to the KPZ universality class ( [Majumdar and Comtet 2004[[Majumdar and Comtet 2005| 
(see Appendix B.4). Importantly, the steady-state average of the extreme advance 
scales with the habitat size in the same fashion as the width w, i.e., (Amax) ~ L" with 
q:r;1/2 [Fig. [^b)]. That is, the expected displacement of the front-runner ahead of the 
invader's mean position increases with habitat size in the above power-law fashion. 
Again, certain system-speciflc characteristics still depend on local rates, neighborhood 
size, etc., but our general ecological prediction does not depend on these details. 

Figure [71^a,b) provides an example, using two different sets of parameters. As 
habitat size increases, the mode and dispersion of the estimated probability densities 
P{Amax,L) increase, driven by the dependence of {Amax) on L. Rescaling these posi- 
tively skewed histograms, and plotting the resulting distributions of {Amax / {Amax}) , 
collapses the data articulately for all system sizes and both sets of parameters onto the 
(theoretical) scaled Airy density [Figure [7Kc,d)]. Note that no fits are used in collapsing 
the scaled data with the theoretical curve. The largest deviations from the theoretical 
prediction can be observed around the peak of the distributions for some larger sys- 
tem sizes. This poorer statistical sampling is, in part, due to the correlated nature of 
the steady-state time series of the extremes, which is progressively more significant for 
larger system sizes. 

Longer fronts increase the expected maximal difference between locations of the 
front runner and the front's mean position. The scaling of the front runner's rel- 
ative position is possible because the distribution of the extremes is governed by 
this single scaling {Amax) ~ L ' . Note that traditional extreme-value limit theo- 
rems (|Gumbel 19281 IFisher 19281 IBerman 1964|) for weakly- or un-correlated fiuctua- 
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Fig. 7 (a) Steady-state distributions of the extreme advance, relative to mean front, for several 
habitat sizes. Parameters are ai = 0.5, 02 = 0.7, and fi = 0.2 ("Set 1"), and (b) for 01 = 0.7, 
02 = 0.8, and /i = 0.1 ("Set 2"). (c) Scaled histograms for all habitat sizes, and for both Set 1 
and Set 2. The solid curve is the scaled Airy density, Eq. 1IB.8II (|Majumdar and Comtet 2004[ |. 
(d) Same data as in (c) but on lin-log scales. 



tions along the front would yield only a weak logarithmic increase with the habitat size 
L ( [Raychaudhuri ct al. 2001|IGuclu and Korniss 2004) . 

To emphasize the significance of universality, we repeated some of our analyses for 
the extremes of the front and compared our model's behavior with with two other, 
well-known invasion models, the Eden model and the basic contact process. We note 
that the fronts in both of these models have long been believed to belong to the 
KPZ class (JuUicn and Botet 1985a: JuUicn and Botct 1985bl lPlischke and Racz 19851 
IKertesz and Wolf 1988 Moro 2001 Ferrcira and Alvcs 2006ll . The scaled probability 
densities exhibit the expected, progressive data collapse on to the universal Airy distri- 
bution [Fig. [8]. The scaled extreme-distribution plots include not only different models 
and system sizes, but for illustration, different neighborhood sizes as well (5 — 4, 8, 12 
for the contact process). The Eden model, the basic contact process, and our model 
of invasion under preemptive competition share scaling properties, revealing common. 




Eclm. L=100 
a Ellen. L=200 
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Fig. 8 (a) Scaled steady-state distributions of the extreme advance for the contact process 
(CP) and for the Eden model for L = 100, 200, 400. For both models there is only one species 
present (invaders) with local propagation rate 02 and mortality rate /x. For the CP, we set 02=1 
and /i=0.20. The front dynamics of the CP for «2=1 and /i=0 is equivalent to a variant of the 
Eden model llJullien an d Botet 1985a! OuUien and Botet 1985bt . and we used these parameter 
values to obtain the Eden data. For the Eden model, we used only nearest-neighbor propagation 
(5=4). For the CP, we also explored different neighborhood sizes, (5 = 4, 8, 12; all arc shown on 
the plot. The solid curve is the scaled Airy density, Eq. I IB.8II jMajumdar and~Comtet 2004| ). 
(b) Same data as in (a) but on lin-log scales. 



key properties of their advancing fronts, in the hght of KPZ universality. Consequently, 
we can infer a series of ecological predictions shared by the most often cited discrete, 
stochastic models of invasive growth. 



5 Discussion 



We recently applied theory for nucleation in homogeneous environments to the spatial 
dynamics of invasion (|Korniss and Caraco 2005l|O'Malley et al. 2005||Q'Malley et al. 2006a[ ). 
Immigration introduces individuals of an initially rare, but competitively superior, in- 
vader randomly independently in space and time. An introduced invader, stochastically, 
may die or may propagate locally and initiate a cluster of invaders. A nucleation event 
occurs when an invader cluster reaches a critical size where further growth through 
propagation and decline through mortality have equal probability (jGandhi et al. 19991 
lYasi et al. 2006IIAllstadt et al. 2007|l . Supercritical clusters, on average, exhibit radial 
symmetry and expand at a constant velocity after sufficiently long periods; their growth 
excludes the resident competitively. Circular fronts will reach the same final velocity as 
their linear counterparts, with the same corrections as times increases (O'Malley et al. 
in press) , so that our analysis of linear fronts furthers application of nucleation theory 
to the evolutionary ecology of invasion. 

To characterize important consequences of discreteness and stochasticity for in- 
vasive advance, we introduced kinetic roughening of the invader-resident interface. 
Roughening of our simulated fronts behaves as a member of the KPZ-universality 
clciss in one (transverse) dimension. Since steady-state roughening depends on habitat 
size, the distribution of the invader's maximal incursion, relative to the front's mean 
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position, also depends on habitat size. The advancing front's steady-state velocity cor- 
rection varies inversely with habitat size, again a consequence of the two-dimensional 
stochastic front's roughening. 

Kawasaki et al. (2006) report two-dimensional stochastic simulations of an invader 
advancing into empty space; propagation into an open site can occur if any of its 8 near- 
est neighbors is occupied. They also assume that invaders are immortal. These assump- 
tions, with growth initiated from a compact region occupied by the invaders, render 
their model essentially a variant of the Eden model (|Jullien and Botet 1985a")l . Hence 
their model should exhibit KPZ roughening (|Jullien and Botet 1985b1IPlischke and Racz 19851 
IKertesz and Wolf 19881 IFerreira a nd Alves 20061) . 

Kawasaki et al. (2006) suggest that stochasticity accelerates invasive advance in 
their two-dimensional simulations, and that reports of deceleration due to stochastic- 
ity may arise because the latter address only one-dimensional processes ( [Snyder 2003[ ). 
Above, we detailed simulation results finding that stochasticity (combined with dis- 
creteness) significantly slows the velocity of an invading competitor in two-dimensions. 
Kawasaki et al.'s (2006) seemingly atypical result reflects the particular deterministic 
model they compared to their simulations, rather than dimensionality. We compared 
our stochastic invasion velocities to results from the corresponding reaction-diffusion 
model, the common baseline approach to invasion velocity in ecology (lAndow et al. 19901 
IHolmes et al. 1994|l . In the context of Moro's (2003) analysis, our stochastic model has 
1 < A'^* < 10; i.e. less than 10 individuals in each row of the interface width. Our 
deterministic comparison, the reaction-diffusion model, has A'^* -^ oo. Our stochastic 
velocity is slower due to diffusion limitation; our pushed fronts are slower than a pulled 
front. Kawasaki et al. (2006) base their claim on comparing stochastic simulations to 
simultaneous advance of the invader in each row. The invader advances (periodically) 
as an exactly straight "cliff;" no roughening can occur. Their deterministic compar- 
ison behaves as if N* — >■ 0, since the front has no width; that is, Kawasaki et al. 
(2006) compare their stochastic simulations to an idealized pushed front, rather than 
to a pulled front (both being inadequate to estimate the velocity of strongly-correlated 
rough fronts). 

Antonovics et al. (2006) present a "patch model" where local demes of a single 
species are coupled by movement between neighboring locations. A gradient in density- 
independent mortality restricts the species' range, and the density profile predicted by 
a mean field model refiects the spatial pattern in mortality. Stochastic simulations 
with small deme sizes produce average density profiles falling below the mean field 
prediction; increased rates of diffusion diminish effects of stochasticity, and density 
profiles approach the mean field result (| Antonovics et al. 2006|l . Hence, a reduction 
in diffusion limitation may accelerate an advancing front or extend a static range 
boundary. 

Our model assumes local propagation only, to simulate limited dispersal. Den- 
sity profiles at the advancing front strongly constrain mixing of the invader and resi- 
dent along the interface. The resulting diffusion limitation (|Moro 200ip distinguishes 
reaction-diffusion velocity from invasion velocities driven strictly by local births, bud- 
ding, or vegetative propagation. The very narrow region with a non-zero density gra- 
dient in any row (along the direction of the invader's advance) of our model's interface 
indicates little mixing of the two species, and compares well with Holway's (1998) study 
of the invasive Argentine ant. New colonies lie close to existing colonies. Budding of 
the invader colonies pushes native ants back as the front advances, and the width of 
the front between invader and residents apparently remains quite narrow. 
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Although our model has no explicit diffusion, the coarse-graining procedure does 
generate, albeit nonlinear, diffusion in the continuum limit (see Appendix A). Including 
further local diffusive moves in our model will not change its universality class, provided 
that the effective diffusion coefficient is not much larger than the effective reaction rate 
HMoro 200ip . Thus, one way to change the asymptotic scaling and the corresponding 
universality class of the front is to introduce local diffusion with large rates (compared 
to local propagation), or alternatively, introduce long-range dispersal. In these cases, 
one expects a crossover from the KPZ class to that of the deterministic reaction- 
diffusion equation, or similar model (|Kot et al. 1996[l . which exhibits a pulled front. 

The dynamical crossover from pushed to pulled fronts may help organize results 
on the modelling and measurement of invasive growth. Generally, we can associate 
different fronts with ecological-invasion studies categorized according to the partic- 
ular spatial scale emphasized. At the within-habitat scale, some plant and animal 
invaders' spatial growth is driven by neighborhood-level propagation (|DAntonio 19931 
ISilvertown et al. 1994l[Schwinning and Parsons 1996||Holway 1998] ). The limitation on 
empirically observed dispersal distances appears consistent with effective-diffusion lim- 
itation and pushed fronts, within any habitat invaded successfully. Jump dispersal, at 
the between-habitat scale, for these same species then must depend on other ecological 
or cultural processes generating introduction events (jKorniss and Caraco 2005]! . 

As estimated dispersal distances of individuals increase (lAndow et al. 19901 |Bjornstad et al. 2002[ ) 
diffusion limitation relaxes, and a reaction-diffusion dynamics may approximate spatial 
advance reasonably well, at both the within-habitat ( |Dwyer 1992|IFrantzen and van den Bosch 2000|l 
and between-habitat (|van den Bosch et al. 19921 INash et al. 1995")) scales. Reaction- 
diffusion models, as noted above, predict pulled fronts advancing at a constant asymp- 
totic velocity. Another set of invasion models can admit increasing velocities; they 
emphasize the between-habitat spatial scale and note the importance of long-distance 
dispersal for some biological invasions IjKot et al. 19961 ICannas et al. 2006;) . Particular 
applications include the biogeographical range expansion of tree species, as indicated by 
paleoecological records l|Clark et al. 1998]) , aerial spread of pathogens at the continen- 
tal scale (| Aylor 2003 I , and dispersal vectored by human migratory (jFisher et al. 200T|) 



or economic (|Ruiz et al. 2000]) activity. Long-distance dispersal events, even if rare, 
can establish colonies detached from the main body of invaders; the two populations 
may then grow and coalesce ( [Shigesada et al. 1995[ ). Invasion dynamics averaging over 
dispersal-distance distributions associated with these models can predict pulled fronts 
with an ever-increasing spatial velocity; Clark et al. (2001) discuss details and theoret- 
ical modifications. 



5.1 The front-runner 

Our front-runner analyses integrate a series of stochastic spatial models through a 
few universal scaling relationships. For all models belonging to the KPZ universality 
class (in one transverse dimension) and for large enough habitat size L, {Amax)^L ' 
[Fig. I2b)]; the distribution of the scaled probability density is given by the Airy dis- 
tribution [Figs. [7]and[S]. These universal features will be exhibited by any invasion 
model belonging to the KPZ class (|Plischke and Racz 19851 IFerreira and Alves 2006p . 
not just our two-species model, the contact process, and the Eden model. 

We referred to transverse system size Ly — L as habitat size, not ordinarily a "tune- 
able parameter" for probing roughness characteristics of the interface. Furthermore, for 
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simplicity, we employed periodic boundary conditions in the transverse direction. From 
the viewpoint of empirical application, we can envision a scenario where one has access 
to only a limited observation window on the full landscape. Designate the window's size 
-^obsi where 1 <C ^obs *^ L. These length scales are measured in units of the minimum 
linear size needed to sustain a single individual (linear size of a lattice site in simula- 
tions). One of the research challenges in extreme- value statistics is to predict some char- 
acteristic of the extremes on global scales (i.e., at the entire habitat or landscape level) 
based on some limited local measurements (i.e., taken from the observation window). 
One can divide the observation region into a number of smaller non-overlapping sub- 
windows of size I, such that 1 <^l <^ -i-obs) smd by changing the size of the subwindow I, 
perform a finite-size analysis to extract estimates for the roughness exponent from the 
average width and the extremes as w{l) ~ Z" and {Amaxipj) ~ I" ■ The subwindows will 
now have "window" or free boundary conditions (|Antal et al. 20011 lAntal et al. 20021 
[Majumdar and Comtet 2005 [ ), clearly differing from the periodic boundaries we simu- 
lated. However, boundary conditions do not affect the roughness exponent, hence the 
scaling properties of width and that of the extremes with I. Therefore, one can estimate 
the size of the extremes for the full landscape of size L from measurements obtained in a 
limited observation window, {Amax[L)) « (Z\rnaa;(I'obs))(i/iobs)"- In contrast, were 
the interfacial fluctuations weakly- or un-correlated, traditional extreme- value statistics 
(for short-tailed parent distributions) would apply. In turn, extrapolating local mea- 
surements to global scales for the extremes would yield only marginal (logarithmic) 
increase, {Amax{L)) « (Zlmaa;(iobs))[ln(-I')/ln(Lobs)]^^''', where 7 is a non-universal 
exponent, and is related the precise shape of the exponential-like tail of the parent 
densities ( |Raychaudhuri et al. 200T]lGuclu and Korniss 2004[l . 

The velocity of spatial advance and the properties of the interface region remain 
objects of study in population biology because of their conceptual significance and 
potential importance to applied ecology. Despite attention directed to invasive spatial 
growth, ecological theory has overlooked the impact of correlations along an advancing 
front (roughening in a two-dimensional environment). Although these spatial correla- 
tions arise through local interactions, at a scale sometimes considered "noise," they 
can directly affect population/community-level phenomena: invader velocity, dynamics 
at the interface, and relative position of the front-runner. 
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Appendix 

A Deterministic reaction-diffusion limit and pulled fronts 

O'Malley et al. (2006b) write the master equation eorresponding to transition rates in Ex- 
pression 1^ (van Karnpen 1981[|. A ssuming that densities at different sites remain uncorre- 
lated ( McKanc and Newman 20041 IKorniss 1997< , one obtains the dynamics of the ensemble- 
averaged local densities pi(x, t)=(ni(x, t)}, 

Pi{x,t + 1) - Pi (x, t) = [1 - pi (x, t) " P2 (x, t)] 

Xy X! P«(x'.*) - A^PiCx,*) , (A.l) 

x'enn(x) 

where i = 1, 2 for resident and invader species, respectively, and we let 5=4 for simplicity. 
Taking the naive continuum limit of the above equations yields the (coarse-grained) equations 
of motion 

^^i^ = ^[1-Pi (X, t) - p2 (X, t)] V V. (X, t) 
at 4 

+ ai[l - pi (x, t) - p2 (x, t)] Pi (x, t) 

- fj.pi{x,t) , (A.2) 

i = 1, 2, and V^ is the Laplacian operator. In statistical physics, equations of this kind are 
referred to as "mean-field" theory, in the sense that correlations between spatial degrees of 
freedom are neglected tMcKane and Newman 2004IIAntonovics et al. 20061 1. This contrasts to 
population biology's terminology where "mean-field" ordinarily refers to non-spatial models. 
After adding appropriate noise terms to the above equations, the resulting Langevin equation 
can, in principle, effectively capture both spatial and stochastic effects, identical to those in 
the underlying discrete individual-based model llSchmittmann and Zia 19951 IHinrichsen 2000l 
Ivan Kampen 1981 1 [Gardiner 1985l l. 

The spatially homogeneous solutions of the above equations, (pj|,p2), are (0,0), (1 — 
/i/«i,0), and (0, 1 — p/a2); coexistence is not feasible. In the parameter regime of interest, 
fj, < ai < «2i only the last solution (0, 1 — p/a2) is stable. The advance of a front separating 
the stable (0, 1 — p-/ct2) (invader dominated) and unstable (1 — p/ai, 0) (resident dominated) 
regions amounts to propagation into an unstable state IIFisher 19371 |Kolmogorov et al. 1937] 
[Aronson and W einberger 1978[ l , a phenomenon that has generated a vast amount of literature 
l]van Saarloos 2 0031 1. Given the mean-field equations IIA.21 1 , the front is "pulled" by the leading 
edge into the unstable state. For a sufficiently sharp initial density profile ( [Murray 2003^ , the 
asymptotic velocity is determined by the infinitesimally small density of invaders that intrude 
into the linearly unstable region dominated by the resident species. Linearizing Eqs. ( IA.2t 
about the unstable state, pi = 1 — p/«i +</>!, P2 = </'2, one obtains for the density of invaders 



/'2(x,f) 

at 



^— V202(x,t)-Hpf— -l')02(x,t). (A.3) 
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Performing standard analysis ( [Murray 200 3 van Sa arloos 200311 of this equation with D = 
(/i/4)(a2/ai) and r = p{a2/ai — 1), the effective diffusion and reaction coefficients, respec- 
tively, we obtain the asymptotic velocity of the "marginally stable" invading fronts, Vp = 
2\/I5r, i.e., 

Vp = — V«2("2 - ai). (A. 4) 

ai 

Vp is the minimal velocity of a traveling wave permitted by Eqs. ( IA.2I I. and is actually realized 
by the deterministic nonlinear reaction-diffusion dynamics for sufficiently sharp initial pro- 
files ([Murray 2003 van Saarloos 20031. Note that we also performed the numerical iteration 
of the nonlinear difference equations Eqs. dA.ll l [the underlying natural discretization of the cor- 
responding reaction-diffusion system Eqs. (IA.2|| ] (O'Malley et al. 2006b O'Mal ley et al. 2009[ |. 
The results show strong agreement with velocity of the linearized equations Eq. (IA.4I I IFig.lSlb)] . 
Despite the nonstandard nonlinear (but deterministic) system of Eqs. (IA.2II . the front velocity 
can be fully governed by the linearized system where the front is pulled by the leading edge. 



20 

In discrete individual-based models with limited diffusion (resulting in strongly corre- 
lated front fluctuations), one does not expect Eq. I IA.41 1 to provide a good approximation 
for the asymptotic front velocity IIMoro 20011 IMoro 2003T l. Instead, Eq. llA.411 can serve as a 
baseline reference for comparison with the actual front velocity. For the deterministic front 
just discussed, where invasion is dominated by infinitesimal densities at the leading edge, 
one can also obtain the width of the interfacial region l |Murray 2003| Ivan Saarloos 200311 . 

w ~ WD/r = (1/2)(1 — ai/a2)~^'^ , independently of the landscape size Ly. 

In passing, we consider the velocity results for our model's reaction-diffusion approximation 
in light of a spatial-competition model by Lewis et al. (2002). Our model assumes preemp- 
tive competition. Since resident and invader use the limiting resource so similarly (interact 
strongly), our model's mean-field approximation prohibits coexistence. The invader's advance 
translates one boundary equilibrium into another. As noted above, iterating Eqs. l lA.ll l pro- 
duces a velocity matching the pulled-front velocity Vp predicted by the linearizing our model's 
reaction-diffusion approximation. That is, for all parameter combinations we evaluated, the 
"linear conjecture" (MoUison a nd Levin 1995l l predicts the mean-field velocity. 

Lewis et al. (2002) write a two-species competition model of Lotka-Volterra form, with each 
species diffusing along a one-dimensional continuum. The spatially homogeneous equilibria are 
standard; coexistence is stable if self-regulation in both species is stronger than interspecific 
competition. The authors derive a set of sufficient conditions under which the linear conjec- 
ture will predict their model's invasion velocity. Outside parameter ranges specified by these 
conditions, their model's velocity might or might not match that obtained by linearizing the 
dynamics at the invasive front's leading edge. Put simply, they find that the linear conjecture 
holds when both diffusion of the competitively superior invader is not too small compared to 
the resident's diffusion, and the species do not interact too strongly. Our results. Fig. |[5b). 
are numerical; if the result for Lotka-Volterra competition with diffusion IILewis et al. 2002l l 
applies more generally, our mean-field approximation's velocity might or might not match 
Vp so closely for other parameter values. We emphasize strongly that our important results 
on invasion velocity concern the differences between the discrete, stochastic model's velocity 
(a pushed front) and that of the associated reaction-diffusion model (a pulled front), quite 
independently of whether or not the latter can be approximated via the linear conjecture. 



B The KPZ equation and scaling behavior of the KPZ universaHty class 

The interface dynamics of a large class of discrete individual-based models can be described 
by the effective (or coarse-grained) Langevin equation IIKardar et al. 1986t 

^^.^ + A(^)%C(..t) (B.l) 

where u is the effective "surface tension," A is the nonlinear coupling constant, and f (j/, t) 
is delta-correlated Gaussian noise, (C(2/i *)) = a^nd (C(s/i t)C(y') *')) = 2F5(j/ — y')S{t — t'), 
with r being the noise intensity. 5{. . .) is the Dirac delta. The coupling constants and the 
noise intensity depend strongly on the details of the course-graining procedure (i.e., the local 
rates or the finite neighborhood size in the original discrete model, and how the continuum 
limit was taken from that model). The variable h{y,t) is the coarse-grained (or mesoscopic) 
limit of the local front advancements hy{t). In some cases, this procedure can be carried out 
explicitly from first principles llPlischke et al. 19871 iKorniss et al. 200^ IKorniss et al. 200311 . 
More often, however, due to the complexity of the details of the model, this task is insur- 
mountable. Even then, however, fundamental symmetry considerations can strongly motivate 
the applicability of the KPZ equation. Thus, equations of this sort should not be interpreted 
as a rigorous continuum limit derived from the local interface rules, but rather as a coarse- 
grained description. This approach has been enormously successful in advancing understanding 
of large-scale morphological properties of interfaces and surfaces in complex natural and arti- 
ficial systems jBarabasi and Stanley 1995||Halpin-Healy and Zhang 1995| IKorniss et al. 2000l 
IBru et al. 200311 . A number of models in the large-scale limit can be effectively described by 
Eq. l IB.lll . i.e., interesting attributes of these models (e.g., the width, the extreme advance, 
and their distributions) exhibit the same scaling behavior and shape as given by the solution 
of Eq. IIB.ll l. These models are said to be belong to the KPZ universality class. 
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B.l Temporal and habitat-size scaling of KPZ fronts 

The temporal and habitat-size scaling of the width (ui^(L,i)) typically identifies the univer- 
sality class of an advancing front. In finite systems, the width grows as (w^{L,t))^t'^l^ from 
early to intermediate times. At a system-size-dependent crossover time, tx~L^, it saturates 
(reaches steady state) and scales as {w^^^) = (w-^(L,oo)) ~ L-^" , where L is the transverse 
linear system size, a, /3, and z are referred to as the roughness, growth, and dynamic expo- 
nents, respectively. The exponents obey the scaling relationship a = /Sz. The full temporal and 
finite-size behavior (reproducing the growth and the steady state regimes as limiting cases) 
can be captured by the Family- Vicsek scaling form | |Family and Vicsek 1985^ 

{w^{L,t))=L^''f{t/Ln. (B.2) 

For small values of its argument, f{x) behaves as a power law, while for large arguments it 

approaches a constant 

f(\—f ^'^^ for X < 1 ,„ „, 

^^ >~ \ const, for a; > 1 ' ^ ■^> 

yielding the scaling behavior of the width, first in the growth regime and then in the steady- 
state regime [Eq. lO, Fig. [3], provided the scaling relationship a = j3z holds. For the KPZ 
universality class, these exponents can be obtained exactly, a = 1/2, j3 = 1/3, and z = 3/2 
lI Kardar et al. 1986 Barabasi and Stanley 19951). Furthermore, in the steady state, employing 
path-integral techniques, both the width distribution (Fol tin et al. 199 41 and the distribution 
of the extremes (Majumdar and Comtct 2004) were obtained analytically. Scaling behaviors 
then can be compared to those obtained from a discrete model (e.g., by simulations), indicating 
whether or not the discrete stochastic model belongs to the KPZ universality class. 



B.2 Steady-state width distribution of the KPZ universality class 

The width distribution typically provides a strong signature of the underlying universality class 
of the fiuctuating, growing interface IIAntal et al. 200T1 lAntal et al. 200211 . It is important to 
note that for models with local dispersal, the measured probability density (or histogram) of 
the individual (row y) interface fluctuations (or "parent" densities) P(hy — h) have tails that 
decay faster than any power law (approximately Gaussian in our model, but this local property 
depends strongly on model details). We also note that if the height fluctuations {hy — h)^ in 
Eq. l[2]l were independent (or their correlation along the interface decayed sufficiently fast), 
the probability density of w^ would be Gaussian, by the central limit theorem. Instead, for our 
roughened front, these variables are strongly correlated. As universality arguments suggest, 
the shape of the steady-state distribution of the width for all (sufficiently large) system sizes 
and for all models belonging to the KPZ class will be identical, governed by a single scale. 
That is, the average width {w^) ~ L. Hence, the width distribution for any model in this class 
with rough fronts can be written as 

P(«;2,L) = -i-<P(«;2/(^2)), (B.4) 

where ^(s) is the probability density of the scaled variable s = w'^/{w^). 'P{s) for periodic 
boundary conditions (used in our simulations) was obtained by Foltin et al. (1994), and is 
given by 

2 °° 
0(s) = — ^(-l)"-in2e-(-'/6)"'» . (B.5) 

The universal scaling function above can be employed to test whether a discrete individual- 
based invasion model with a rough front belongs to the KPZ class. Analyzing the width his- 
tograms, as empirical measures of the theoretical probability density, supports our conclusion 
that the two-species invasion model, the focus of our investigations, belongs to the KPZ class 
[Fig. [9]. Note that no fits were used in collapsing the scaled data to the theoretical curve. 
Analytic scaling functions for the width distribution can also be obtained for other types of 
boundary conditions, e.g., "window" or free IIAntal et al. 20021 . 
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Fig. 9 (a) Steady-state width distributions in our two-species invasive front for several habitat 
sizes with parameters oi=0.5, 02=0.7, /i=0.2 ("Set 1"); and (b) for ai=0.7, 02=0.8, ^=0.1 
("Set 2"). (c) Scaled histograms for all habitat sizes, for both Set 1 and Set 2. The solid curve 
is the scaled analytic width distribution of the KPZ class, Eq. IIB.St lIFoltin et al. 1994ll . (d) 
Same scaled data as in (c), but on lin-log scales. 



To emphasize the concept of universality, we also analyzed data and constructed his- 
tograms for propagating fronts in two well-known invasion models, the Eden model and the 
contact process. We note that the fronts in both of these individual-based models have long 
been believed to belong to the KPZ class (JuUicn and Botct 1985a JuUicn and Bot et 1985bl 
IPlischke and Racz 1985IIKertesz and Wolf 1988 Fcrrcira and Alves 2006 Moro 2001i) . The scaled 
probability densities exhibit progressive data collapse onto the universal scaled KPZ density 
[Fig. 1101 , offering further support that the invasive fronts in the Eden model and in the basic 
contact process belong to the KPZ class. 



B.3 Asymptotic velocity of KPZ fronts: temporal and habitat-size corrections 



As described in the text, a model's asymptotic front velocity cannot be universal. However, 
Krug and Meakin (1990) showed that the forms of the temporal and finite-size {Ly = L) 
corrections are universal within a given class. Above we demonstrated that the propagating 
interface in our two-species invasion model belongs to the KPZ universality class with ex- 
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Fig. 10 (a) Scaled steady-state width distributions for the contact process (CP) and for the 
Eden model for L = 100, 200, 400. For both models there is only one species present (invaders) 
with local propagation rate «2 and mortality rate fi. For the CP, we set Q2=l and ^=0.20. 
The front dynamics of the CP for 02 =1 and fi=0 is equivalent to a variant of the Eden 
model t JuUien and Botet 1985al l. and we used these parameter values to obtain the Eden 
data. For the Eden model, we used only nearest- neighbor propagation (5=4). For the CP, we 
also explored different neighborhood sizes, (5 = 4, 8, 12; all are shown on the plot. The solid 
curve is the scaled analytic KPZ width distribution, Eq. IIB.SI l IIFoltin et al. 1994l l. (b) Same 
scaled data as in (a) but on lin-log scales. 



poncnts a = 1/2, /3 = 1/3, and z = 3/2. Using Krug and Meakin's result, one obtains the 
corrections to the front velocity v{t, L), 



v{t,L) : 



* 


- cit-2(i-Qr)/z = „. _ cit-2/3 for t < L^ 


* 


- C2L-2{l-<«) =v* - C2L-1 for t > L^ 



(B.6) 



where ci and C2 are non-universal constants that depend, as does v* , on propagation and 
mortality rates at the level of individual competitors. 

The first expression in Eqs. IIB.6II indicates that the early temporal correction to the front's 
velocity scales as 0{t~^'^). Pulled fronts produced in deterministic reaction-diffusion models 
exhibit 0{t~^) corrections to the asymptotic velocity IIAndow et al. 19901 IHolmes et al. 19941 1. 
A discrete model's pushed front has a larger temporal correction; effects of time since initiation 
of invasion (from t=0 until tx^L^) persist longer in the stochastic case. In the simulations 
we measure h{t)/t as it follows the same temporal and system-size scaling as those of v{t, L) 
[Eq. I IB.6I 1I (O' Malley et al. 2006b^ . Figure [TTT a) supports the above scaling for our discrete 
individual-based two-species invasion front. 

More importantly, the second expression in Eq. 1IB.6II predicts that after roughening equi- 
librates, the invading species' velocity must be corrected according to the size of the habitat 
invaded. This correction scales as 0{L~^). Steady-state velocity increases with habitat size; 
longer invading fronts move faster. Figure [TTT b) and its inset show how this asymptotic be- 
havior is controlled by system size L. Observed corrections to the velocity v* decline linearly 
as a function of L~^, in accordance with the second expression Eq. 1IB.6II . The physical the- 
ory for roughened interfaces explains the effect of habitat size on frontal velocity, and further 
gives the form of the quantitative dependence. These corrections to the steady-state velocity 
were observed qualitatively for the basic contact process (Ellnc r~et al. 1998ll and for the Eden 
model UKawasaki et al. 200611 . but neither study addressed either the quantitative correction 
for habitat (system) size or the underlying basis for the dynamic behavior: KPZ roughening 
of the front. 
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Fig. 11 (a) Temporal corrections to the asymptotic front velocity for L=800 for different 
values of 02, with oi=0.50 and /i=0.20. The horizontal (time) axis is scaled as t~^i^ , moti- 
vated by the form of the corrections of the one-dimensional KPZ class [Eq. 1IB.6I1 ]. The inset 
enlarges the region of (a) where the universal temporal corrections follow the KPZ behavior. 
The straight solid lines correspond to the linear scaling as a function of t"'^'''^ . (b) Finite-size 
corrections to the asymptotic velocity as a function of L~^ , motivated by the universal cor- 
rections of the one-dimensional KPZ class, for ai=0.50, O2=0.70, and /i=0.20. The straight 
dashed line corresponds to the linear behavior as a function of 1/L. 



B.4 Distribution of extremes in steady-state KPZ fronts: the Airy density 

The distribution of the front's extreme advance (relative to the mean) provides an additional 
strong signature of the underlying universality class of a stochastic interface. This distribution 
was only recently obtained in analytic form ( Majumdar and Comtet 2004 , Majumdar and Comtet 2005 | l. 

As universality arguments imply (Schchr and Majumdar 2006), the shape of the steady- 
state distribution of the front's extreme advance (relative to the mean) for all (sufficiently 
large) system sizes and for all models belonging to the KPZ class will be identical, governed 
by a single scale, {Amax) ~ L^'^. For a given model with a rough front one has 



i^yAmax: L/J ■ 



1 



{■^max) 



I'iA^ax/iAmax)) , 



(B.7) 



where !?(«) is the probability density of the scaled variable u = Amax/{Amax)- For peri- 
odic boundary conditions (employed in our simulations) Majumdar and Comtet (2004, 2005) 
obtained the analytic result 

n^) = ^J\fA (yi") , (B.8) 

where fj\ is the Airy density function. A continuous random variable X {X >0) has an Airy 
probability density if: 



fA{X) 



2^6 Y^ 

Xio/3 2^' 



'bUX\2/ijj^^_ 



5/6,4/3,fefc/x2) 



(B.9) 



where fe^ = 2afe''/27, and (7(-, •, •) is the confluent hypergeometric function, a^. is the magnitude 
of the fc-th root of the Airy function Ai{-) on the negative real axis l |Abramowitz and Stegun 1972^ . 
The Airy distribution arises in a remarkable number of different applications, including the 
extremum of a set of correlated random variables; see Majumdar and Comtet (2005), Schehr 
and Majumdar (2006), and Guclu et al. (2007). 
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